function [V_new,cy1,cy2,w1,w2,w_thr,wthr_approx] = Value_Iteration(cy1_old,cy2_old,w1_old,w2_old)

global wgrid betta delta gamma pi1 pi2 theta e1 e2 ey1 ey2 v1aut v2aut wfb vfb nw wmin wmax cy1fb cy2fb

crit=1e-10;
optionsset=optimoptions('fmincon','Display','off','TolX',10^(-8),'TolFun',crit,'Algorithm','sqp');
options=optimoptions('fsolve','Display','off');

%Preallocation Vector
V_new=zeros(nw,1);
cy1=zeros(nw,1);
cy2=zeros(nw,1);
w1=zeros(nw,1);
w2=zeros(nw,1);

%Policy Valuation
PI=zeros(nw); 
    for k=1:nw
    PI(k,:)=PI(k,:)+weight(wgrid,w1_old(k),pi1);
    PI(k,:)=PI(k,:)+weight(wgrid,w2_old(k),pi2);
    end
    
    num=pi1*((betta/delta)*utility(e1-cy1_old,gamma)+utility(cy1_old,gamma))+pi2*((betta/delta)*utility(e2-cy2_old,gamma)+utility(cy2_old,gamma));
    V=(eye(nw)-delta*PI)\num; %Use the value as an educated initial guess
    
   distanza=1;
   
   while distanza>crit 
  Vnew=num+delta*(pi1*interp1(wgrid,V,w1_old,'spline',V(end))+pi2*interp1(wgrid,V,w2_old,'spline',V(end)));
  distanza=max(abs(Vnew-V));
  V=Vnew;
   end
     
%Policy Improvement
  function OBJ_neg = value_unconstrained(c) %c(1): cys1, c(2): cys2, c(3): omegas1, c(4): omegas2 
    
            OBJ_neg=...
            -(pi1.*(utility(c(1),gamma)+(betta./delta).*utility(e1-c(1),gamma)) ...
            + pi2.*(utility(c(2),gamma)+(betta./delta).*utility(e2-c(2),gamma)) + ...
            delta.*(pi1*interp1(wgrid,V,c(3),'spline') + pi2.*interp1(wgrid,V,c(4),'spline')));         
  end

 function [ineqcon,eqcon] = constraints_fun(c,omega)
        eqcon=zeros(1,1);
        ineqcon=zeros(2,1);
        
        eqcon(1)=v2aut-utility(c(2),gamma) - betta*c(4); % Participation constraint s=2 (binding)
        ineqcon(1)=v1aut-utility(c(1),gamma) - betta*c(3); % Participation constraint s=1 (not necessarily binding)
        ineqcon(2)=omega-(pi1*utility(e1-c(1),gamma)+pi2*utility(e2-c(2),gamma)); % Promise keeping constraint      
 end

wthr_approx=max(wgrid(abs(V(1)-V)<1e-8)); 
w_thr=fsolve(@(x)interp1(wgrid,cy1_old,x,'spline')-cy1fb,wthr_approx,options);

lb = [0;0;w_thr;w_thr];
ub = [ey1;ey2;wmax;wmax];  
startguess=[cy1_old cy2_old w1_old w2_old];

for k=1:nw
    guess=startguess(k,:)';
    [out,V_new(k,1)]=fmincon(@(x)value_unconstrained(x),guess,[],[],[],[],lb,ub,@(x)constraints_fun(x,wgrid(k)),optionsset);
    cy1(k,1)=out(1);
    cy2(k,1)=out(2);
    w1(k,1)=out(3);
    w2(k,1)=out(4);
    
end
  V_new=-V_new;
end

